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ABSTRACT 

Context. Circumstellar disks are considered to be the environment for the formation of planets. The growth of dust grains in these 
disks is the first step in the core accretion-gas capture planet formation scenario. Indicators and evidence of disk evolution can be 
traced in spatially resolved images and the spectral energy distribution (SED) of these objects. 

Aims. We develop a model for the dust phase of the edge-on oriented circumstellar disk of the Butterfly Star which allows one to fit 
observed multi-wavelength images and the SED simultaneously. 

Methods. Our model is based on spatially resolved high angular resolution observations at 1.3 mm, 894 [im, 2.07 [im, 1.87 [im, 
1.60 [im, and 1.13 [im and an extensively covered SED ranging from 12 [im to 2.7 mm, including a detailed spectrum obtained with the 
Spitzer Space Telescope in the range from 12 [im to 38 [im. A parameter study based on a grid search method involving the detailed 
analysis of every parameter was performed to constrain the disk parameters and find the best-fit model for the independent observa- 
tions. The individual observations were modeled simultaneously, using our continuum radiative transfer code. 
Results. We derived a model that is capable of reproducing all of the observations of the disk at the same time. We find quantitative 
evidence for grain growth up to ~ 100 [im- sized particles, vertical settling of larger dust grains toward the disk midplane, and radial 
segregation of the latter toward the central star. Within our best-fit model the large grains have a distribution with a scale height of 
3.7 AU at 100 AU and a radial extent of 175 AU compared to a hydrostatic scale height of 6.9 AU at 100 AU and an outer disk radius 
of 300 AU. Our results are consistent with current theoretical models for the evolution of circumstellar disks and the early stages of 
planet formation. 

Key words, protoplanetary disks - stars: pre-main sequence - stars: individual: IRAS 04302+2247 - circumstellar matter - planets 
and satellites: formation - radiative transfer 



1. Introduction 

Circumstellar disks are a natural outcome of the star forma- 
tion process and are known to dissipate over time (~10 6 -10 7 yr; 
Hai sch et al.|2001| ) by means of stellar winds or by photoevapo- 
ration caused by either the radia tion of the central star or by an 
external so urce of radiation (e.g.,|Hollenbach et al.|2000| [Clarke 
|et al.|200T] [Alexander & A rmitage 2007), by accretion onto the 
central star ([Hartman n et al.||1998|), by grain growth and frag- 
menta tion (e.g., |D ullemond & Dominik 2005; Dominik et al. 



2007 



2011 



Natta et al.||2007; 



Garaud et al. 2012 



tion of planets (e.g., Pol 
|2007] |2&T0l [Fortier et al. 



Birnstiel et al.fefTJ (Sauter & Wolf 



Ubach et al. 2012), and by the forma 



ack et al.||1996[ |Boss 2002; Armitage 
2012| ). The way in which the disk ma- 



terial dissipates has important implications for the possibility of 
planet formation. 

During the last decades several theories to explain the for- 
mation of planets in circumstellar disks have been pro posed, 
for example, the core accretion-gas capture scenario ([Pollack 



et al. 1996; Papaloizou & Terquem 2006 ; |Lissauer & Stevenson 



2007). Sub-micron- sized dust particles, which are strongly cou- 
pled to the motion of the gas in the disk via gas drag, agglomer- 
ate through low- velocity collisions to eventually form planetes- 
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imals of several kilometers in size (e. g., |Beckwith et al. 2000; 
Dominik et al.|2007[ [Natta et al.|2007| ). During this coagulation 
process, larger particles decouple from the turbulent gas motion, 
settle toward the disk midplane and radially drift toward the in- 
ner parts of the disk as a result of the impact of stellar gravity 
and gas drag (e.g., Weidenschilling 1977 ; Barriere-Fouchet et al. 
2005l|Fromang & PapaloizoupOO^HD'Alessio et al.|2006| ). The 
growth beyond pianetesimais occurs via direct collisions, with 
an increasing role for gravitational focusing as masses become 
larger, leading to the gravitat ional agglomer ation of these bod- 
ies to rocky planets (Safronov & Zvjagina 1969). A phase of 
runaway growth occurs, followed by an oligarchic growth phase 
dKokubo & Ida| 19981 |Thommes et al.|2003] ). 

There are, however, many unresolved problems in the picture 
of this formation process. For example, it is uncertain how dus t 
grains overcome the radial drift barrier ( [Weiden schilling 1 1977 ). 
This occurs as it is expected that meter- sized bodies migrate to- 
ward the star in a timescale that is shorter than the timescale 
for further growth which effectively stops further grain growth. 
Another example is the fragmentation barrier, where the boul- 
ders are destroyed at typical collisi on speeds, halting the dust 
growth at centimeter to meter sizes (Benz 2000 ; Blum & Wurm 
2008 J |Brauer et al.||2008| ). Moreover, the dust coagulation pro- 
cess, the first stage of planet formation, must also overcome the 
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charging b arrier (Okuzumi 2009 ) and bouncing barrier ( |Zsom 
|et al.|2010| ) which affect the later stages, too. 

Besides the hypothesis that planetesimals may form from 
pairwise collisional growth of smaller bodies alone, there is the 
hypothesis that, when the dust particles have reached cm- scale 
as a consequence of this collisional growth, the planetesimals 
may form from the gravitational fragmentation of a dense par- 
ticle sub-disk near the equatorial plane resulting from the dust 
settling. Protoplanetary disks may support regions within which 
turbulence acts to locally enhance the ratio of solids to gas. As a 
result, patches of the disk may become dense enough to trig- 
ger a gravitational instability and collapse into planetesimals 
( Safronov & Zvjagina 1969; Goldreich & Ward 1973; Armitage 



2007; Chiang & Youdin 2010| ). This suggestion entirely by- 
passes the size scales that are most vulnerable to radial drift and 
would allow the radial drift and fragmentation barrier to be over- 
come. 

With the current observational capacities it is not possible 
to directly observe large dust boulders (> 1 m) and planetesi- 
mals, but signs and implications of their formation can be de- 
tected. The stratified structure resulting from grain growth and 
settling has a significant imp act on observable quantities of the 
disk (e.g., Dullemond & Dominik 2004). The spectral energy 
distribution (SED) provides constraints on the dust mass and 
grain-size distribution. Scattered light images at near-infrared 
(NIR) and mid-infrared (MIR) wavelengths probe the surface 
layers of the optically thick disks and reveal properties of small 
dust grain s through wavelength-dependent opacity (e.g., Watson] 
|et al.|2007] ). Continuum observations at submillimeter and mil- 
limeter (mm) wavelengths are mainly sensitive to the bulk of the 
disk mass closer to the midplane and play an important role in 
constraining disk structure proper ties, for example, the spatial 
dust distribution of larger grains ( Andrews & Williams 2005 , 
2007 ). The analysis of scattered light or re-emission images or 
the SED only provides a limited view of a disk leading to ambi- 



guities in the derived disk and dust properties (e.g., Chiang et al. 
2001 ). To precisely study the physical processes like dust evo 



lution and to strongly reduce model degeneracies, it is essential 
to combine resolved images at different wavelengths and a SED 
with a good wavelength coverage in a multi- wavelength obser- 
vational and modeling approach. 

In this paper we study the circumstellar disk of the promi- 
nent Butterfly Star, also known as IRAS 04302+2247, usin g this 
approach. This Class I young stellar object (YSO; e.g., | Adams] 
|et al.||1987[ |Lada 1987) is located in the Taurus- Auriga molec- 
ular c loud complex at a distance of ~ 140 pc (Ken yon et al. 
1994] ). It is surrounded by an edge-on seen flared circumstel- 



lar disk, so the centr al star is not visib l e directly (inclin ation 
i = 90° ± 3°; |Padgett et"aT][T999l [200T| |Wolf et al.|2003) . As 



shown by interferometric observations at submm and mm wave- 
lengths, the disk exhibits a radius of about 300 AU ([Wolf et aT 



2003 2008| ). Several attempts have been undertaken to model 



the structure and physical conditions in the circumstellar envi- 
ronment of this object (e.g., ILucas & Roche|| 19971 |Wolf et al 



2003 Stark et al.|[2006). In particular, Wolf et al. (2003) were 



the first to base their modeling on scattered light images and re- 
solved mm maps and, therefore, were the first to consider wave- 
length regimes tracing different physical processes in different 
regions of the circumstellar environment. They concluded that 
the grains in the outer parts of the circumstellar environment are 
comparable to those from the interstellar medium (ISM), while 
dust grains have grown via coagulation by several orders of mag- 
nitude in the much denser parts of the disk. However, the obser- 
vational data they presented did not allow them to constrain the 



spatial dependence of the dust grain properties in the disk inte- 
rior. 

We present new observations and a coherent multi- 
wavelength model for the circumstellar disk of the Butterfly Star 
that accounts for spatially resolved data sets ranging from the 
NIR to mm wavelengths as well as for the well- sampled SED of 
this object. These observations provide different, complemen- 
tary detailed perspectives on the disk structure and dust proper- 
ties. Here, coherent means that the model is capable of repro- 
ducing all considered observations obtained at different wave- 
lengths at the same time. This new modeling became necessary 
as new observations are not properly reproduced by the model 
of |Wolf et al.| ( |2003l |2008| ). Furthermore, these new observa- 
tions with much higher angular resolution promise to allow one 
to put stronger constraints on the disk structure and dust grain 
properties. 

The observations and data reduction are introduced in 
Sect.[2] In Sect.[3| we characterize the data set that is the basis 
for our modeling. A detailed description of our model can be 
found in Sect.]?] and the entire modeling process is described in 
Sect.[5] The results of this study are presented in Sects.[6]and|7 
and their implications are discussed in Sect.[S] Finally, Sect " 
contains a summary and concluding remarks. 



2. Observations and data reduction 

The data set used in this study includes both spatially resolved 
images at mm, submm, and NIR wavelengths, as well as a well- 
sampled SED. In this section these observations and the data 
reduction procedures are discussed. 



2.1. Millimeter observation 

The Butterfly Star was observed with the IRAM Plateau de Bure 
interferometer (PdBI, [Lucas| 1 99 1 1 ) in its two most extended con- 
figurations (A and B). The dual polarization receivers were tuned 
near 220 GHz (A = 1364 jim). The wideband WIDEX correlator 
provided a total bandwidth of 4 GHz in each of the two polariza- 
tions. The covered frequency range included lines of the 13 CO 
and C 18 isotopologues, as well as H 2 CO and SO. Spectral line 
results will be reported in a separate paper, but the continuum 
data discussed here was generated by avoiding line contamina- 
tion. 

Observations were performed on February 3 and 4, 2010 
(Aq configuration, two baselines of 730-760 m, excellent RMS 
(baseline-based) phase noise < 35°), January 22 (Aq configu- 
ration), and February 10, 2011 (B configuration, up to 450 m 
baselines, RMS phase noise < 45°). Although the repeated con- 
figurations produced similar results, data from January 22, 2011 
was not included in the final analysis because of its higher phase 
noise. With very good observing conditions (< 1 mm water va- 
por, going up to 1.5 mm only for the last 2 hours of February 10, 
201 1), single sideband system temperatures ranged between 100 
to 140 K. Phase noise yields an effective seeing of approximately 
0.15". 

Phase calibration was performed using the nearby quasars 
0400+258 and 0507+179, and flux calibration is referred to 
MWC349, with an assumed flux density of 1.7 Jy at this fre- 
quency, with an overall calibration accuracy of order 10 %. 

With natural weighting, the angular resolution of the data is 
0.59" x 0.35" at a position angle (PA) of 31.1°. The shortest 
baseline is large, 79 m including projection effects, so that struc- 
tures larger than about 4" could be substantially resolved out. 
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However, given the source elongation and orientation and the 
beam shape, deconvolution is efficient in recovering most struc- 
tures. 

The total flux density is about 101 mJy ± lOmJy. The ex- 
pected thermal noise in the map is 0. 1 8 mJy/beam, but because 
of phase noise, the dynamic range is limited and the effective 
noise level in the deconvolved map is ~ 0.4 mJy/beam. 

2.2. Submillimeter observation 

Observations of the Butterfly Star with the Submillimeter Array 
(SMA, |Ho et al.|2004| ) were carried out on January 9, 2006 us- 
ing the upper and lower sidebands at 330 GHz and 340 GHz, re- 
spectively (A = 894 |im). The angular resolution of the data is 
0.67" x 0.53" at PA = -76.5°. The SMA observations and the 
data reduction were described in detail by |Wolf et al.| p008 ). 

2.3. Near-infrared observations 

The NIR observations of the Butterfly Star were obtained 
with the Near Infrared Camera and Multi-Object Spectrometer 
(NICMOS, |Thompson et~aL] [T998] ) on the Hubble Space 
Telescope (HST) on August 19, 1997. The NICMOS NIR 
Camera, NIC2, and the filters F110W (A c = 1.13 |im), F160W 
(A c = 1.60 urn), F187W (A c = 1.87 urn), and F205W (A c = 
2.07 |im) were used. The data reduction and flux calibration were 
described in detail by Padge tt et al.| ( |1999| ). Figure [3] shows a 
composite image of the Butterfly Star at NIR wavelengths. 



1.3 mm map shows a clear maximum at the center, a local mini- 
mum is visible at this position in the 894 |im map. The nature of 
this local minimum is discussed in Sect. 18.11 




Fig. 1. 1.3 mm map obtained with the PdBI. Negative contours 
are dashed and at -3 <x, where cr is the background noise of the 
map. The positive contour levels are in steps of 5 cr starting at 3 cr 
and going up to 53 cr. The ellipse in the bottom left illustrates the 
shape and orientation of the beam. 



2.4. Spectroscopic data 

The Butterfly Star was observed with the low-resolution 
modules of the Spitzer Space Telescopes (SST) InfraRed 
Spectrograph (IRS) on October 1, 2007 and March 1, 2004 for 
short-low (SL) and long-low (LL), respectively. The LL obser- 
vation was obtained as part of the IRS guaranteed time obser- 
vations (PI: Houck), while the deeper SL observations were ob- 
tained as part of the open time program 30765 (PI: Stapelfeldt). 
The total exposure time was 587 seconds for the SL modules and 
12.6 seconds for the LL modules. The two-dimensional droop- 
corrected frames obtained using pipeline version SI 6. 1.0 were 
co-added for each nod position, and the two nods were pairwise 
differenced to remove the background. One-dimensional spectra 
were extracted using apertures of 3 pixels for the SL modules, 
while 4 and 5 pixel apertures were used for the LL2 and LL1 
modules, respectively. The one-dimensional spectra were flux- 
calibrated using the spectral response function included in the 
SST pipeline. 




3. Basis for modeling 

The basis of our modeling forms the (sub)mm and NIR maps, as 
well as the SED, including the S ST/IRS spectrum. This data set 
is briefly characterized in this section. 



Fig. 2. 894 |im map obtained with the SMA. The contour levels 
are in steps of 2 cr starting at -4 cr, leaving out the zero contour 
level, and going up to 16 cr. The ellipse in the bottom left illus- 
trates the shape and orientation of the beam. 



3.1. Images 

The interferometric continuum maps at 1.3 mm and 894 |im are 
shown in Figs. [T] and [2] The emission seen in these maps is ther- 
mal re-emission of the dust. While the source is unresolved ver- 
tically to the disk in both maps, it is well-resolved in radial direc- 
tion. The circumstellar disk seen in both maps has a north/south 
extension of ~ 600 AU. While the brightness structure of the 



The emission seen in the NIR images is scattered light from 
the central star (see Fig.[3]below and Fig. 1 in Wolf et al. (2003)). 
Their appearance is dominated by a totally opaque linear feature 
that bisects the scattered light structure. This feature is inter- 
preted as a circumstellar disk seen edge-on. Because of the large 
optical depth of this disk, no point source is detected at any of 
the observed NIR wavelengths. The dependence of the height 
of the disk on the wavelength is clearly visible. The bipolar 
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scattered light nebula, which extends ~ 900 AU in the direction 
north/south, also shows a complex morphology with approxi- 
mately equal brightness between the eastern and western lobes. 
As this appearance reminds one of a butterfly, the name Butterfly 
Star has been established ( [Lucas & Roche|1997] ). Comparing the 
(sub)mm maps with the NIR images it can be seen that the posi- 
tion and orientation of the elongated (sub)mm structure perfectly 
fits that of the dark lane in the NIR. 

We rotate the observed images by the position angle of the 
image z/-axis and the position angle of the Butterfly Star (PA = 
175 °, as simulated images are symmetric we use PA = -5 °) in 
order to align the major axis of the dust lane with the vertical 
axis, as it is for the simulated images. 



Table 1. Photometric data points for the Butterfly Star used in 
the modeling 




Fig. 3. Composite image of the Butterfly Star seen in scattered 



light. For details see Sects. 3.1 and 7.2 



3.2. Spectral energy distribution 

The photometric measurements used in this modeling are pre- 
sented in Table[T] In addition, the S ST/IRS spectrum in the range 
from 12 |im to 38 |im is considered. In our disk modeling we 
do not include the scattered-light part of the SED below 12 [im 
as it is not an objective of this study to reproduce t he h ighly- 
structured envelope seen in the NIR images (see Sect.|5.2|). 



Figure [4] shows the complete SED of the Butterfly Star con- 
sisting of the data shown in Table[T] the S ST/IRS spectrum rang- 
ing from 5 jim to 38 \im, as well as HST/NICMOS, SST/IRAC, 
and WISE data points. The SED of our target is typical for a 
Class I YSO. It peaks close to 100 [im and shows a 10 \im ab- 
sorption feature. This silicate feature and the steep slope of the 
IRS spectrum beyond 10 \im require a high inclination angle 
which is consistent with the spatially resolved NIR images. In 
Fig. [5] the fit to the (sub)mm slope of the SED of the Butterfly 
Star can be seen which yields the millimeter spectral index 
<*mm = -dlog(F /l )/dlog(/l) = 2.48 ± 0.07. This value is signif- 
icantly smaller than that found for small ISM-sized dust grains 
(-3.7), indicating grain growth in the disk (see, e.g., Natta et al. 
12007] ). 
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Fig. 4. Spectral energy distribution of the Butterfly Star. 

4. Model description 

4.1. The disk 

Our model for the Butterfly Star consists of a parameterized disk. 
We describe this disk with a radially and vertically dependent 
dens ity distribution based on the work of |Shakura & Sunyaev| 
( |1973| ) which can be written as 



Aiisk =Po\—J ex Pl"2 



h(r cy i) 



(1) 



Here, z is the cylindrical coordinate with z = corresponding to 
the disk midplane and r cy i is the radial distance fro m this z-axis 
(see, e.g., [Wolf et al.||2053} |Stapelfeldt eFaLl^998l see Fig. |6] 
for illustration). The parameter p is determined by the total dust 
mass mdust of the disk. The quantity ro is a reference radius with 
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I observation 
fit 



cv = 2.48 ±0.07 



400 
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Fig. 5. Spectral index a mm of the Butterfly Star. For details see 
SectEl 



ro = 100 AU, whereas h is the vertical scale height and a function 

of r cy i 



h(r cy i) = h 



(2) 



The parameters a and J3 describe the radial density profile and 
the disk flaring, respectively. The disk extends from an inner 
cylindrical radius r m to an outer one r out . Together with the 
quantity ho, these five parameters are used to adjust the disk 
structure in order to fit the data. This ansatz has already been 
successfully applied to other edge-on circumstellar disks, such 
as HH30 (Mad lener et"aTl|20l2) , CB 26 jSauter et al.1 2009), 
IMLupi (Pinte et al.|2008])7HV Tauri C ( [Stapelfeldt et~aL|2 003), 
and HK Tauri ( [Stapelfeldt et aE1[T998] ). Integrating Eq. [1] along 
the z-axis yields the surface density 



£(r cyl ) = S (^ 



with p - a - /5. 



(3) 




Fig. 6. Illustration of a flared edge-on oriented disk (xz plane) 
described by Eq. [T] The contour levels are logarithmic (log 10 ) 
with outwards decreasing values. 



4.2. Dust 

In our model we focus on radiative transfer through the dust 
that domina tes the transport of radiation and thermal structure 
of the disk ( Chiang & Goldreich |1997| ). The density structure 
described in Eq. [l]is given by the gas in the disk. Because the 
dust in the disk is coupled to the gas, its density distribution is 
also described by this equation. The dust grain properties in our 
model can be divided into three groups: the shape of the dust 
grains, their chemical composition, and their size distribution. 



4.2.1. Grain shape 

We assume the dust grains to be homogeneous, spherical, non- 
aligned, and non-oriented particles (Mie theory), although they 
are expected to feature a much more co mplex and fractal struc- 
ture. As shown by Voshchinnikov (2002), shape, chemical com- 
position, and size of dust grains cannot be determined separately, 
but only in combination. Hence, we restrict our model to the sim- 
plest but least ambiguous shape possible. 



4.2.2. Chemical composition 

To model the chemical composition of the dust grains we use 
the homogeneous mixture of smoothed astronomical silicate and 
graphite with a mean de nsity of p gra in = 2.5 g cm" 3 and optical 
properties described by Weingartner & Draine| ( |2001| ). Relative 
abundances of 62.5% astronomical silicate and 37.5% graphite 
are used ( [Draine & Lee|1984[|Weingartner & Draine|200l] ). This 
grain model h as been employe d successfull y in the mo deling of 
HH 30 ( |Madlener et"aTJ|20ni ) and CB 26 ( |Sauter et al.]|20U9l ), 
for example. Given the index of refraction, the optical attributes 
of a homogeneous sphere can be calculated using Mie theory. 
We extrapolate the complex refractive indices to a wavelength 
of 2.7 mm to cover all available observational data. This is ap- 
plicable since both the real and the imaginary part of the refrac- 
tive index show asymptotic behavior in this wavelength regime. 
Because graphite is a highly anisotropic material, it is necessary 
to treat the two possible alignments of the electric field to the 
crystal axis independently using the so-called | - | approxima- 
tion for graphite spheres. That means, if gext is the extinction 
coefficient, then 



1 2 

Gext = ^Gext(^ll) + ^GextOJ 



(4) 



explains the dependence of gext on the components of the dielec- 
tric tensor (e\\ and e ± ) of the electri c field parallel and perpen- 
dicular to the crystallographic axis. Draine & Malhotra] ( 1993 ) 
showed that this approximation is sufficiently accurate for ex- 
tinction curve modeling. 

4.2.3. Grain sizes 

The sizes of the dust grains in our model are distributed accord- 
ing to a power law of the form 



dn (a) - a 3 5 da with 



< a < a n 



(5) 



Here, a is the dust grain radius and n(a) the number of dust 
grains with a specific radius. For = 5 nm and a max = 250 nm 
this grain- size distribution becomes the commonly-known MRN 
distribution of the ISM by |Mathis et aT] ( |1977| ). To properly 
model the different grain sizes in a dust grain mixture, an ar- 
bitrary number of separate dust grain si zes w i thin a g iven in- 
terval [amin : <2 max ] has to be considered. |Wolf| ( |2003a| ) showed 
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that the observables resulting from radiative transfer (RT) sim- 
ulations considering each grain species separately are close to 
those based on weighted mean dust grain parameters of the dust 
grain mixture. Thus, we use weighted mean values for the op- 
tical properties of the dust grain mixture. Moreover, since the 
separate processes (e.g., grain growth, dust settling, grain-grain 
interaction, mixing processes) and their mutual influence dur- 
ing the evolution of the circumstellar environment are still rather 
poorly understood and also to reduce the number of free param- 
eters necessary to describe the dust grain mixture, we assume no 
spatial dependence of its properties and assume the power law 
distribution (Eq.[5]) to be valid throughout the circumstellar disk. 

4.3. Heating sources 

Stellar heating is the primary heating source for the circumstellar 
disk. The radiation of the central star heats the dust which then, 
in turn, re-emits at longer wavelengths. The star, which is rep- 
resented by a single black body, is characterized in our model 
by its effective temperature and bolometric luminosity L*. 
Owing to obscuration by the circumstellar disk, the illuminat- 
ing source of the system cannot be observed directly, leading to 
weak constraints on these two parameters from observation. 

Viscous heating of the disk is neglected and we ignore accre- 
tion and turbulent processes in the disk. 



5. Means of modeling 

The goal of the modeling process is to find a coherent model 
for the circumstellar disk of the Butterfly Star explaining all the 
described observational data. Using continuum RT simulations 
we compute observables from the model that are then compared 
to the observed data to find the best-fit model. In this section we 
describe the means of modeling in detail. 

5.1. Radiative transfer 

For our continuum RT simulations we use the program MC3D 
dWolf et al.|1999[|Wolf|2Q03bl ). It is based on the Monte-Carlo 
method and solves the continuum radiative transfer problem self- 
consistently. It makes use of the temperature correction tech- 
nique described by Bjork man & Wood] ( |2001| ), the absorption 
concept introduced by Lucy| (T999|), and the enforc ed scatter- 
ing scheme proposed by |Cashwell & Everett ( 1959). The opti- 
cal properties of the dust grains (scattering, extinction and ab- 
sorption cross sections, scattering phase function) and their in- 
teraction with the radiation field is calculated using Mie theory. 
Multiple and anisotropic scattering is considered. 

In order to derive a spatially resolved dust temperature dis- 
tribution, the model space has to be subdivided into volume el- 
ements inside which a constant temperature is assumed. Both 
the symmetry of the density distribution and the density gradient 
distribution have to be taken into account. We use a spherical 
model space, centered on the illuminating star and an equidis- 
tant subdivision of the model in ^-direction, while a logarithmic 
radial subdivision is applied to resolve the temperature gradient 
at the very dense inner region of the disk. 

The RT is simulated at 107 wavelengths; 100 of them are 
logarithmically distributed in the wavelength range [0.05 |im, 
2000 [Am] and the remaining 7 wavelengths are distributed in the 
range [894 [im, 3000 |im]. The quantities we derive with MC3D 
from the model are 



2. NIR images at 1.13 |im, 1.60 |im, 1.87 |im, and 2.07 |im; and 

3. TheSED. 



5.2. Model fitting properties 

The simulated (sub)mm maps are convolved with the point 
spread function (PSF) that is described by an elliptical Gaussian 
function based on the beam major and minor axis as well as the 
beam position angle from the corresponding observation (see 
Sect. [2]). To account for the rotation of the observed maps, the 
beam position angle of the PSF is adjusted by the same angle. 
For the simulated NIR images we use the corresponding PSF 
obtained with the PSF modeling tool Tiny Tim v.7.4 ( |Krist et aL] 
201 l[|Krist & Hook|2004] ) for the convolution. Additionally, the 
convolved simulated maps have the same pixel scale as the ob- 
served data. This way we ensure that the simulations are compa- 
rable to the observations. 

There are primary features that we want to reproduce with 
our model. These also determine the criteria for the best-fit 
model. For the SED we aim to reproduce the thermal re-emission 
spectrum over almost three orders of magnitude from the MIR 
down to the mm wavelength regime. For resolved images the is- 
sue is not as simple as for the SED. Our model is rotationally 
symmetric and thus does not account for any related asymmetry 
that is seen in the observations. We consider the maps of the disk 
to be two different groups, the (sub)mm maps and the NIR maps. 



• Maps in the (sub)mm: Although the two (sub)mm maps are 
simply structured they provide two decisive features that 
constrain our model and that we want to reproduce, the peak 
flux density and the spatial brightness distribution (mainly 
the radial brightness distribution with constraints for the ver- 
tical extent). The uv data is, in general, used to analyze the 
(sub)mm data because this avoids the non-linear deconvolu- 
tion step. Here, as the source is quite strong and barely re- 
solved in one direction, the deconvolution is straightforward 
and we use image plane comparison. 

• Maps in the NIR: The four images in the NIR show more 
structures and details than the (sub)mm maps. In addition 
to the disk, which appears as a dark lane, a very complex 
and highly-structured envelope surrounds the disk. It is not 
an objective of this study to reproduce this envelope with its 
wavelength-dependent morphology, so we restrict ourselves 
to reproducing the width of the dust lane and, consequently, 
the wavelength dependence of the width of the dust lane. 
This is the main information provided by the NIR images be- 
cause it constrains the vertical opacity structure of the disk. 

5.3. Quality of the fit 

For each comparison between model and observation on the 
above points, we determine a goodness of fit value ^ that char- 
acterizes how well the model fit matches the observational data 
and so the best-fit model is represented by the smallest goodness 
of fit value of the investigated parameter space. 

SED: The goodness of fit of the SED ^ ED is given by 



^SED 



1 N 

-Y 



err 



(6) 



(Sub)mm maps at 894 \xm and 1.3 mm; 



Here, N is the number of observed wavelengths, cot are the ob- 
served flux densities and fa the synthetic flux densities at these 
wavelengths. The observational uncertainties are considered by 
<Ti. As the re-emission part of the SED of the Butterfly Star, 
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which we want to fit, is better sampled in the wavelength range 
from 12 |im to 38 |im thanks to the S ST/IRS spectrum than at 
longer wavelengths up to 2.7 mm, we include weighting factors 
gi in the calculation of the £s ED . They are choosen in such a 
way to ensure that the better and the worse sampled parts of the 
SED are equally considered in the £ 2 ED calculation, i.e., that the 
flux measurements ranging from 12 |im to 38 |im have in total 
the same weighting as the data points in the range from 70 |im 
to 2730 |im altogether. A natural weighting fit would rely ex- 
cessively on a narrow range of wavelengths and could possibly 
miss gross features which have a much lower S/N. By using an 
adjusted weighting, we ensure that we do not sacrifice the overall 
agreement for tiny details in the SED. 

Maps: For every map the goodness of fit £^ ap is calculated 

by 



2 _ (Jd-0J) 
£map - 2 



(7) 



where \i is the modeled data, co the observed data, and <x the ob- 
servational uncertainty. Therefore, for the (sub)mm maps every 
pixel j whose flux density in the observed data is larger than the 
background noise in a square box with a side length according 
to the diameter of the object and centered on the object is taken 
into account, which lets us rewrite Eq.[7]as 



(sub)mm map 



(8) 



In this case, TV is the number of pixels taken into account in this 
box and <x represents the background noise of the observed map. 
Instead of a map-based goodness of fit there is also the possibil- 
ity of choosing an approach where only the profile along the disk 
axis is taken into account. However, this method would not allow 
for any direct constraints on the disk scale height. Therefore, we 
decided to choose the approach described above 

To get a total goodness of fit %\ 
fit we combine all ^ 



t 2 otal that characterizes a model 



£ 2 

S total 



1 N 

k= 1 



(9) 



Based on the £ 2 otal we get from Eq. [5J we give our modeling 
errors, i.e., constraints on the model parameters, as the range 
where we can alter the parameter values without changing £ 2 otal 
by more than 10%. Allowing for a larger variation of £ 2 otal than 
10% gives generally worse results. 



5.4. Parameter space study 

The modeling process is divided into two steps. We simultane- 
ously fit the (sub)mm data for the first time and figure out if it 
is possible to find a model fit that reproduces these data using 
the described model (Sect.|4]) with a single dust grain-size dis- 
tribution. Based on the results of this first modeling step (see 
Sect. |6.2| ) we modify our model setup and also include the SED 
and the NIR images in the fitting process. 

On the basis of the model described in the previous sec- 
tions, the parameter space we deal with is defined by several 
free parameters. These adjustable model parameters are used to 
repr oduc e the key features of our observations as described in 
Sect.|5.2| The volume of the parameter space consists of six and 



nine free parameters in the first and second modeling step, re- 
spectively. For our study, the selection and range of these param- 
eters is based on modeling of other similar objects and previou s 
modeling efforts (see, e.g., |Sauter et al.|200 9; Wolf et al.|2003] ). 
To constrain the model parameters as strongly as possible and to 
find the parameter set of the best-fit model we use a grid search 
based method. 



6. First modeling step 

We first consider a disk model with no spatial variations of the 
dust grain properties. 



6.1. Model Parameters 

The six adjustable model parameters we use in the first step are 
the following: 



The exponents a and ft that describe the radial density profile 
of the disk and the disk flaring, respectively (see Eqs.[T]and 
[2]). In our modeling, both parameters are treated as indepen- 
dent quantities. We vary a in the range of [1.0, 4.0] and /3 in 
the range of [1.0,2.0]. 

The scaling factor ho of the disk's scale height (see Eq. [2]). 
Because the value of the reference radius ro is 100 AU, ho 
equals the scale height at a radial distance from the star of 
r cy i = 100 AU (h = fc(100 AU)). We consider h in the range 
of [1.0, 25.0] AU. 

The dust mass md us t of the disk. The dust mass range we 
probe is [1.0 x 10" 5 , 1.0 x 10" 2 ] M Q . 

The inner disk radius r m . Inner holes in protoplanetary disks 
are a common prediction of theories of disk evolution and 
are observed in various circumstellar disks (e.g., Alexander 
& Armitage||2009l |Andrews et al.|[20TT| |Grafe et al.||2011 



Esp aillat et al.|2012| ). The inner radius is varied in the range 



of [0.1, 50.0] AU. The lower limit is approximately the dust 
sublimation radius. 
• The maximum grain radius a mSLX . Grain growth is a major 
issue for protoplanetary disks as it is the first step toward the 
formation of planets. For this upper radius of the grain- size 
distribution we consider nine different values: [0.25, 1.0, 10, 
30, 50, 75, 100, 300, 1000] jim. 

For the minimum grain radius a m { n we use 5 nm for both mod- 
eling steps. The outer radius of the disk r out is fixed at 300 AU 
for both steps bas ed on the (sub)mm observations and the mod- 
eling results from |Wolf et al.|p003| ). On the basis of the mass 
of the star of ~ 1 .7 M Q resulting from observations of the gas 
kinematics (Dutrey e t al., in prep.) and pre-main - sequence evo- 
lutionary tracks (e.g., Hillenbrand & White 2004), we found that 
L* = 5 L© and T+ = 4500 K characterize the central star best. 
These values are compatible with the radiative transfer simula- 
tions and are used in the whole modeling. According to the pre- 
cisely determined inclination i of the circumstellar disk of the 
Butterfly Star of 90° ± 3° (Sect.[j}, we fix the inclination at 90° 
in the first modeling step. Because of its location in the Taurus- 
Auriga molecular cloud complex, a distance d of the object of 
140 pc is adopted throughout the modeling process. An illustra- 
tion of the density structure of a model fit from step one is seen 
in Fig. [6] Table|2] shows an overview of the fixed parameters used 
in the modeling process. 
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6.2. Results 

The results presented in this subsection are based on the model 
and applied assumptions outlined in Sects.[4]and[5] 

Using the described modeling setup of the first step, which 
is mainly characterized by a single dust grain-size distribution, it 
is not possible to find a coherent model that properly reproduces 
the (sub)mm data. Figures [7] and 8] show the best-fit results for a 
grain-size distribution with small ISM-sized dust grains (a max = 
0.25 |im) and larger dust grains (a max = 100 |im), that yield the 
motivation for the second modeling step. Dust grains with sizes 
comparable to those found in the ISM (a max - 0.25 \xm) are 
neither capable of reproducing the brightness structure of the 
(sub)mm maps nor do they allow us to fit the s hallo w mm slope 
of the SED (Fig. [7] see also Fig. [5] and Sect. |3.2| ). The radial 
brightness distributions of the (sub)mm data in Fig. [7] show an 
insufficient flux density at almost all regions in the disk. Only 
in the outermost disk regions the small dust grains yield a flux 
density that is in agreement with the observation. From the SED 
in Fig. [7] it can be seen that the MIR and far-infrared (FIR) part 
is reproduced using ISM-sized dust grains, but that it is not pos- 
sible to account for the (sub)mm part of the SED. These results 
strongly suggest the need for larger dust grains in the disk, par- 
ticularly in the region close to the central star. 

We find that the brightness structure of the (sub)mm maps is 
reproduced best using a dust grain- size distribution with ci mdLX — 
100 |im (Fig. [8]). Furthermore, the large grains fit the (sub)mm 
part of the SED better than the small ones while particles signif- 
icantly smaller or larger than a max = 100 \xm are inapplicable. In 
contrast, it is not possible to account for the IR part of the SED 
using the large dust grains. Moreover, they yield a flux density 
that is too high, especially in the outer disk regions as can be 
seen in the (sub)mm radial brightness profiles (Fig. [8]). This in- 
dicates that small dust grains, especially in the outer parts of the 
disk, are still needed in the model. 

Nevertheless, with no combination of the used variable pa- 
rameters it is possible to account for the proper shape of the ra- 
dial brightness profiles. We conclude, that to properly fit all the 
observational data and to find a coherent model that explains all 
these data, we get from the first modeling step that we need large 
dust grains with a max -100 \xm as well as small ISM-sized dust 
grains at the same time in the disk model. Additionally, from the 
disk scale height /z(100 AU) that amounts to 6 AU and 12 AU us- 
ing the small and the large dust grains, respectively, and from 
the described brightness profiles (Figs. [7] and [8]), we conclude 
that we need the large grains to be located around the disk mid- 
plane with a small vertical extent and concentrated toward the 
central star. 

In addition, no need for a hole in the inner disk is found from 
the modeling. As a consequence we fixed the inner disk radius 
in the second modeling step at 0.1 AU which is approximately 
the dust sublimation radius. 



7. Second modeling step 

The results of the first modeling step show that it is not possible 
to properly fit the (sub)mm data at the same time and that at least 
two dust grain-size distributions are necessary, with larger dust 
grains concentrated around the disk midplane. 



Table 2. Overview of the fixed parameters. 



Parameter 


Value 


a min [nm] 


5 


r in [AU] 


0.1 


rout [AU] 


300 


U [L Q ] 


5 


[K] 


4500 


d [pc] 


140 



Notes. The inner disk radius r in is only fixed in the second modeling 
step. 



adjustable parameters in the second step instead of six. In con- 
trast to the first step, we consider two different dust grain- size 
distributions in the disk model. Their location in the disk can 
be seen in Fig. [9] They differ in the maximum grain size with 
^max,id > ^max,sd (Id: large dust, sd: small dust). The parameter 
^max,sd is set to 0.25 |im which corresponds to the du st grain-size 
distribution found in the ISM flMathis et aL]|1977| ) and that we 
expect in the outer regions and upper layers of the circumstel- 
lar disk of the Butterfly Star ( |Wolf et aL|2 003). Because it is not 
necessary to vary the inner disk radius r m we fix it at 0. 1 AU. The 
minimum grain radius a^ n for both grain-size distributions, the 
outer radius of the disk r out , the distance d of the object, L*, and 
T+ are the same as in the first modeling step (see also Table [2]). 
The parameters a, ft, ho, and md us t are used again as free param- 
eters. In addition, we use the following adjustable parameters in 
the second step of our modeling: 



• The two quantities fa and r ou t,id that describe the vertical and 
radial extent of the region where the grain- size distribution 
with the larger dust grains is located (see Fig. [9]). The param- 
eter fa is radially dependent in the same way as in Eq. [2] 
We emphasize that fa is not a scale height in the sense of 
ho. In our model the circumstellar disk is described by one 
density distribution (see Eq.[T]). Therefore, the vertical extent 
described by fa has a sharp upper limit. The range we con- 
sider for fa (at 100 AU) is [1.0, 25.0] AU and for r ouUd it is 
[50, 300] AU. 

• A scaling factor sm that adjusts the distribution of the dust 
mass of the disk within the two different dust regions. This 
parameter is affected by the other free disk parameters, i.e., if 
the total number of dust grains is changed in either of the two 
different dust regions by any of the adjustable parameters, 
the ratio of the dust mass of these two regions is modified 
although this scaling factor is not varied. 

• The inclination i of the circumstellar disk. Despite the small 
uncertainty of the inclination measurement, here we allow 
the inclination to vary in the range of 90° ±3°. 

• The maximum grain radius of the dust grain- size distribution 
containing the larger dust grains <z max ,id- We consider three 
different values: [50, 100, 150] jim. 

Using this model setup in the second modeling step we are ac- 
counting for grain evolution and its dependence on the radial and 
vertical position in the circumstellar disk. 



7.1. Model parameters 

Owing to the modification of our mo del setup based on the re- 
sults of the first modeling step (Sect.|6.2|), we use a set of nine 



7.2. Results 

The results presented in this subsection are based on the model 
and applied assumptions outlined in Sects.|4]and[5] 
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Table 3. Overview of the parameter ranges, best-fit values, and constraints on the model parameters of the second modeling step. 



Parameter 


Minimum value 


Maximum value 


Best-fit value 


Lower limit 


Upper Limit 


s 


(X 


1.0 


4.0 


1.2 


1.2 


1.6 


0.2 


6 


1.0 


2.0 


1.14 


1.12 


1.16 


0.01 


h [AU] 


1.0 


25.0 


10.0 


9.0 


11.0 


1.0 


& [AU] 


1.0 


25.0 


6.0 


4.0 


7.0 


1.0 


rouud [AU] 


50 


300 


175 


170 


210 


5.0 


sm 


0.001 


0.9 


0.083 


0.068 a 


0.093° 


o.oor 


^dust [M Q ] 


0.00001 


0.01 


0.0009 


0.0008 


0.001 


0.0001 


a maxM IN 


50 


150 


100 


> 50 


< 150 


50 




87 


93 


90 


90 


90 


1 



Notes. In the last column the step width 6 between the given lower and upper limits of each parameter is given. 

The constraints on the scaling factor sm refer to the best- fit m odel. As this parameter is affected by the other parameters, no constraints for the 
whole investigated parameter space are given (see also Sect.|7.1|>. 





1 out, Id " 



Fig. 9. Illustration of a flared edge-on oriented disk (xz plane) 
using two different grain- size distributions. For details see 
Sect.l7ll 



With the model setup of the second modeling step, which 
is mainly characterized by two different dust grain-size distri- 
butions, it is possible to find a model fit that explains all the 
observational data taken into account and therefore a coherent 
model for the circumstellar disk of the Butterfly Star. The values 
of the parameters of our best-fit model and the constraints on the 
model parameters resulting from our multi-wavelength model- 
ing can be found in Table [3] Only in the small parameter space 
described by the lower and upper limits of the listed parame- 
ters is it possible to find model fits that satisfactorily reproduce 
all observational data, i.e., £ 2 otal does not change by more than 
10 % compared to the best-fit value. Figures [10| and [TT] show 
the radial brightness profiles of the (sub)mm observations over- 
layed with the corresponding counterparts based on the best-fit 
model. The optical depth r at 894 [im amounts to 1.75 and at 
1.3 mm to 0.64. This means that the disk is optically thick at 
submm wavelength and optically thin at mm wavelength. The 



effective specific dust opacity is determined to be 2.58 cm 2 g 1 
and 0.96 cm 2 g" 1 at 894 |im and 1.3 mm, respectively (based on 
the dust model outlined in Sect. |4.2| ). The observed and modeled 
re-emission SED is represented in Fig. 12 In Table [4] the width 
of the dust lane at the point of minimum separation between the 
two lobes in both the observed and modeled scattered light maps 
is listed (see also Fig. [3]). The width of the dust lane stays almost 
constant across the whole radius of the disk such that small devi- 
ations are covered by the listed uncertainties. Overall, our best- 
fit model is in very good agreement with all considered obser- 
vations, especially taking into account the range of wavelengths 
and the variety of observations analyzed in the fitting procedure. 



In Fig. [T0]it can be seen that our model almost perfectly fits 
the 1.3 mm observation. The overall accuracy, i.e., how good the 
modeled profile matches the observed one in terms of the ob- 
served background noise, is 1.39 cr mm + 0.98 cr^. Although the 



central brightness minimum in the submm observation (Fig. 1 1 
is not reproduced by our model we do see a clear flattening of the 
profile w hich is definitely caused by an optical depth effect (see 
also Sect. |8.1| ). With an accuracy of 0.96 <x su b m m ± 0-77 <x su bmm, 
our model also provides a satisfying fit to the submm observa- 
tion. Furthermore, the modeled SED is in very good agreement 
with observations from the MIR over the FIR to the mm regime, 
meaning that almost all modeled flux densities are within the un- 
certainties of the observed flux densities (Fig. [12]). In particular, 
it reproduces the (sub)mm flux densities and hence the mm spec- 
tral index. The width of the dust lane in the modeled NIR images 
is in agreement with the observed width taking into account the 
measurement uncertainties (Table [4}. That means that the four 
modeled NIR images nicely reproduce the observed width of the 
dust lane and therefore also its wavelength dependence which 
was the goal for this wavelength range. 



Table 4. Width of the dust lane in the observed NIR maps and 
those of our best-fit model. 





1.12 [im 


1.60 [im 


1.87 [im 


2.07 [im 


A 


Model 


1.125" 


0.750" 


0.600" 


0.525" 


+0.075" 


Observation 


1.200" 


0.825" 


0.600" 


0.525" 


+0.075" 



Notes. In the last column the uncertainty A of the measurements is 
given. 



10 



C. Grafe et al.: Vertical settling and radial segregation of large dust grains in the circumstellar disk of the Butterfly Star 
0.020 1 1 1 1 1 1 



0.015 



£r o.oio 



10 
£Z 

X 
Z3 



0.005 



0.00Q 



PdBI 

A = 1364 |jm 

a ma x,id = 100 urn 
™ dust = 0.0009 M 
Hq = 10.0 AU 
r outJd = 175.0 AU 
h ld = 6.0 AU 



i - 

a 



= 0.083 
90.0 
= 1.2 



(3= 1.14 



observation 
]— [ with lcr 

errorbars 
- - model 



-"300 -200 -100 100 

Radial distance from the center [AU] 



200 




300 



Fig. 10. Radial brightness profile at 1.3 mm based on the observation and our best-fit model. For details see Sect. 7.2 
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Fig. 11. Radial brightness profile at 894 um based on the observation and our best-fit model. For details see Sect. 7.2 
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Fig. 12. SED based on the observations and our best-fit model. For details see Sect. 
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8. Discussion 



8.1. Brightness minimum in the SMA observation 



The central brightness min imum seen in the observed 894 \xm 

( 2008 ) where two expla- 



map was first described in Wolf et al. 



nations for this minimum are discussed, an inner disk hole and 
an optical depth effect. When comparing the beam size of the 
two (sub)mm maps (see Sect.[2]) it can be seen that the beam 
size of the 1.3 mm observation is significantly smaller than that 
of the 894 |im map. This means that any disk hole that would 
be capable of producing such a local brightness minimum in the 
submm observation would also show up in the 1.3 mm observa- 
tion. However, at 1.3 mm no sign of a central brightness min- 
imum or even a flattening of the profile was found. Therefore, 
this rules out the possibility that the local minimum seen in the 
894 |im map is caused by the lack of emitting dust, i.e., an in- 
ner disk hole. Moreover, when comparing both maps and taking 
the beam sizes into account, we can conclude that the disk of 
the Butterfly Star is optically thick even at submm wavelength. 
Thus, the central brightness minimum in the submm observa- 
tion is produced by an optical depth effect (for details see Wolf 
|et al.|2008] ). This is al so s upported by the optical depth at both 
wavelengths (see Sect. |7.2| ) resulting from our best-fit model. We 
emphasize that the observed brightness minimum has a low sig- 
nificance of ~ 1.4 cr and ~ 1.8 <x with respect to the maxima to 
the left and right of the center, respectively. As our best-fit model 
shows a clearly flattened profile in the submm that accounts for 
the optical depth effec t and is in very good agreement with the 
observation (Sect.|7.2[), it is a satisfying result. 



8.2. Dust evolution 

We found a coherent model capable of explaining all consid- 
ered observations. While it was not possible to achieve this goal 
with a model using a single dust grain-size distribution we suc- 
ceeded by using two different distributions. Both distributions 
differ in the size of the maximum particle radius a max and in the 
spatial distribution in the disk. We found that for a max values of 
0.25 |im and 100 |im for the smaller and larger dust grain- size 
distribution, respectively, fit the data best. The grain size we de- 
termine for the larger dust in the disk is almost three orders of 
magnitude larger than upper grain sizes given for the ISM in 
the literature. This strongly indicates that grain growth is taking 
place in the circumstellar disk of the Butterfly Star. However, 
this dust grain size is still one order of magnitude smaller than 
the maximum grain size derived from modeling other young cir- 
cumstellar disks, such as for IMLupi (Pi nte etal.|20 08 ). There, a 
maximum size of a few mm has been found. In the framework of 
our model we exclude the strong influence of such large grains 
on the (sub)mm observations of the disk of the Butterfly Star. 
Our model clearly shows the need for small ISM- sized as well 
as larger (a max ~ 100 |im) dust grains in this disk. In particu- 
lar, the smaller grains are located in the upper layers and outer 
region of the disk and the larger ones in the inner disk region. 
The findings on the grain sizes are o nly v alid in the context of 
the assumed dust properties (see Sect. A2_ ). Assuming a different 
distribution of the sizes of the dust grains may lead to different 
results. But what in general can be concluded is that the larger 
dust grains have a different spatial distribution than the smaller 
ones. 

Comparing the disk scale height of our model and the corre- 
sponding vertical extent of the region where the large dust grains 
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are located at 100 AU yields that the large dust grains are settled 
toward the disk midplane. Furthermore, comparing the radial ex- 
tent of the large dust region with the outer disk radius it can be 
seen that the large dust grains are not distributed over the whole 
radial extend of the circumstellar disk, but are concentrated to- 
ward the star. 

We conclude that besides grain growth in the disk of the 
Butterfly Star, we find evidence for vertical settling and radial 
segregation of the dust in this disk. These three mechanisms are 
the key issues in the process of planetary core formation, which 
start the formation process of planetesimals. Our coherent model 
not only gives strong constraints on the dust properties in this 
disk, but we can also provide quantitative constraints on the ver- 
tical and radial distribution of the phase with large dust grains 
for the first time. 



of ~ 10 AU at a reference radius of 100 AU is required. The de- 
termined flaring index ft of ~ 1.14 and scale height are compa- 
rable to the values for the YSO IMLupi by |Pinte et aL] ( |2008 ). 
Sau ter et aL] ( |2009| ) also found for the circumstellar disk in the 
Bok globule CB 26 a scale height of ~ 10 AU and with p ~ 1.4 a 
slightly larger flaring index. For the circumstellar disk of HH 30, 



Madlene r et al.| ( [20T2] ) found a flaring index comparable to what 



we find for the disk of the Butterfly Star, but a larger scale height 
(h ~ 15 AU). Our best-fit model has a midplane temperature 
of 20 K at 100 AU. Figure [T3| shows the calculated midplane 
temperature as a function of the radial distance from the cen- 
ter. The determin ed value for T^pp i s very si milar to other YSO s 
such as IMLupi ( [Pinte et al. 12008), CB 26 (|Sauter et al. 2009), 
DGTau, DGTau-b, H LTau (Guilloteau et al.|2011) , and HH30 
flMadlener et al.|2012) . 



8.3. Disk mass 

Our best-fit model shows that we need a dust mass of M^ ust = 
9x 10~ 4 M Q to account for the observed flux densities. This value 
is well constrained (see Table [3]) and stronger variations of the 
dust mass would have a significant influence on the re-emission 
SED, for example. The dust mass of the region where the large 
grains are distributed is 2.2 x 10" 4 M and that of the small grains 
is 6.8 x 10" 4 M . Assuming a typical dust-to-gas ratio of 



Ma 



1 

Too' 



(10) 



we get a total disk mass of ~ 9 x 10 2 M© and determine a 
maximum absorption opacity of 1.97 cm 2 g" 1 at 1.3 mm. This 



is in good agreement with the results of Wolf et al. (2003) and 
very similar to other YSOs su ch as IMLupi (Pinte et al.|2008 ), 
DL Tau, GO Tau, and HL Tau ( [Andrews & Williams |2007] )~ 

The derived disk mass depends on the adopted assump tion s 
about the chemistry and shape of the dust grains (see Sect. |4.2| ). 
Fractal grain shapes and the spatial orientation of every dust 
grain are not taken into account as this would only introduce 
further free parameters without providing a qualitative improve- 
ment of our understanding of the disk properties and distribution 
of large grains. The complex shapes of the dust grains would 
have implications on the light scattering and absorption behav- 
ior of the grains (e.g., |Wright|[T987l |Lumme & Rahola||1994) . 
Besides that, it also can be thought of dust grains with almost 
the same absorption cross section as spherical dust grains but 
with much less mass. Such fluffy particles with a porosity up to 
90% are discussed by Voshchinni kov et al. ( |2007 ). 

Furthermore, we point out that the disk mass is proportional 
to the grain density p gra in and the number of dust grains. Within 
our study we can constrain the disk structure, the grain size, and 
their number, but not the density of one grain. For our investi- 
gation of the circumstellar disk, we used p gra in = 2.5 g cm -3 , but 
a more porous, fractal structure of the dust grains may result in 
a smaller value. In turn, this will alter our estimate for the disk 
mass by the same factor. However, our general results on the dif- 
ferent spatial distributions of the smaller and larger dust grains 
are not affected by these assumptions. 

8.4. Disk structure 



Our multi-wavelength modeling allows us to quantitatively con- 
strain most of the geometrical parameters of the circumstellar 
disk of the Butterfly Star. A flared geometry with a scale height 



10' 
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Fig. 13. Midplane temperature of our best-fit model. Tioo is the 
midplane temperature at 100 AU. 



The hydrostatic scale height H is derived by 



H = 



k B T(r)r 3 
I GM+/jm p ' 



(11) 



_gra\ 

li = 2.33 gmol 1 ( Ruden & Pollack||l99lfr is the mean molec- 
ular weight, and ra p is the proton mass. This yields a value of 
H = 6.9 AU at a radius of 100 AU. To compare this to the ver- 
tical extent of the region where the large dust grains are located 
(fa) an equivalent scale height at a radius of 100 AU can be 
determined by 



h\d = 



' [2 r^ d 

Vtt Jo 



exp|-- 



dz 



(12) 



For our best-fit model this results in a scale height for the large 
dust of h\a = 3.7 AU and for the constraints on ho and fa in a 
range of 2.5 AU< h\& < 4.3 AU. This is substantially smaller 
than fa and substantially smaller than the hydrostatic scale 
height H which strongly indicates that the large dust grains in 
the circumstellar disk of the Butterfly Star are vertically settled. 
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8.4.1 . Surface density 



In Fig. 14 the derived surface density of our best-fit model for the 
circumstellar disk of the Butterfly Star as well as the upper and 
lower limit can be seen. The surface density exponent (see Eq.[5]) 
is found to be in the range of [0.04, 0.48]. These values are lower 
than those predic ted by most theoretical models of disks (p ~ 1 ; 
Bell et aL||1997 ). Furthermore, studies of large samples of cir- 



cumstellar disks in the Taurus-Auriga and Ophiuchus-Scorpius 
star formation regio ns b y, e.g.,[K itamura ^t al.|(|2002] ), |Andrews| 
|& Williams] ( [2007] ), and |Guilloteau et al.| { |2011) show that val- 
ues of p less than ~ 1 are common, but they are in general larger 
than 0.5. Therefore, our result for p is slightly smaller than that 
found for most other objects. We determine the surface density 
to be in the range of [3.3, 15.2] gem -2 and [2.9, 3.6] gem -2 at 
5 AU and 100 AU, respectively, which is consistent with the pre- 
viously mentioned surveys. 



£ of our best-fit model 
upper and lower limit 




10° 10 1 
Radial distance from the center [AU] 



Fig. 14. Surface density. The upper and lower limit results from 
the constraints on the model parameters (see Table [3]). 



8.4.2. Disk stability 

To estimate if the circumstellar disk of the Butterfly Star is self- 
gravitating, we determined the ratio yu sg of the Kepler time tk 
and the local free-fall time Tff 



D" 1 



cyl 



Tff : 




with p c = po 



r cy i 



tk 

Tff 



(13) 



(14) 



(15) 



Here, Qk is the Keplerian angular velocity, G is the gravitational 
constant, and M* is the stellar mass. If p sg <?c 1 the disk is 
dominated by the gravitational potential of the central star and 
therefore Keplerian, i.e., it is non- self-gravitating. In contrast, if 



yu sg » 1 the disk is dominated by the gravitational potential of 
the local mass distribution and so becomes near-Keplerian self- 
gravitating to fully self-gravitating with increasing yu sg . FigurefT5 



shows the radial dependent ratio p sg and it can be seen that in our 
model the disk of the Butterfly Star is non-self-gravitating at all 
radii. This also implies that gravitational instabilities in the disk 
are very unlikely. 
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Fig. 15. Ratio p sg . The upper and lower limit results from the 
constraints on the model parameters (see Table[3]). The transition 
border marks the point where the circumstellar disk becomes 
self-gravitating. 



To check this, we made use of the Toomre criterion ( |Toomre| 
1964 ) that describes when a disk becomes unstable through grav- 
itational collapse. The Toomre Q parameter is defined as 



Cs^K 

ttGi: 



(16) 



where c s is the local sound speed. In terms of Q, a disk is un- 
stable to its own self-gravity if Q < 1, and stable if Q > 1. As 
Fig. [T6] clearly shows, in our model Q > 1 throughout the entire 
disk, which means that the circumstellar disk of the Butterfly 
Star is gravitationally stable at all radii as is expected by the de- 
termination of p sg . 



8.5. Comparison with theoretical models 

We compared our results qualitatively with theoretical models 
for the evolution of dust in circumstellar disks by|Garaud et al.| 
(2004), Garau d & Lin| ( |200^ , |Ba"rriere-Fouchet et aLH2005) , and 
|Gara ud (2007|, for example. The Butterfly Star is characterized 
as a Class I YSO and is surrounded by a circumstellar disk that 
is optically thick even at submm wavelengths and by a substan- 
tial envelope. This implies that this object is rather young and 
potential evolution of dust is in its early phase. The main points 
of current models on the first phase of planet formation are sum- 
marized below. The gas in the protoplanetary disk is partially 
pressure supported, meaning that both the centrifugal force and 
the gas pressure counteract the gravity. Because of an outward 
decreasing gas pressure gradient, the gas orbits at sub-Keplerian 
velocity around the star. Small dust particles are well-coupled 
to the gas. As they do not experience the same radial pressure 
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Fig. 16. Toomre parameter Q. The upper and lower limit results 
from the constraints on the model parameters (see Table [3}. 



gradient as the gas they orbit with the Keplerian velocity and 
therefore feel a net inward force causing them to drift inward. 
As a consequence of processes such as Brownian motion and 
turbulence in the disk, the dust grains agglomerate as a result 
of collisions and form aggregates of ever increasing size and 
mass. In the context of our model we found evidence for the 
early stage of this grain growth phase in the circumstellar disk 
of the Butterfly Star with grain radii of up to ~ 100 \xm. This 
is significantly larger than what is commonly found in the ISM 
(~ 0.25 |im). With the increase of the importance of gravity re- 
sulting from the continuing grain growth, the particles decouple 
from the pure gas motion and settle toward the disk midplane. 
Because they have a faster velocity than the gas does, the larger 
particles see a headwind that saps their angular momentum and 
causes them to spiral in toward the star. This radial drift occurs 
on a time scale that is much shorter than the disk lifetime. From 
our multi-wavelength modeling of the disk of the Butterfly Star 
we also found quantitative evidence for the settling of larger dust 
grains toward the disk midplane while smaller particles are re- 
tained in the upper layers of the disk. We found the larger dust 
grains to be concentrated toward the star which is in agreement 
with the predicted inward spiraling. Moreover, theoretical mod- 
els predict that the gas drag efficiency varies according to the 
grain size, where intermediate- sized particles (100 [im to 10 cm) 
experience the strongest perturbation to their movement. The 
largest particles that we have found based on our model exhibit a 
radius of ~ 100 |im. Therefore, our findings for radial segregation 
of the dust is not unexpected. Theoretical work by [Garaud et al. 



( 2012 ) suggests the coexistence of two particle populations, one 



with larger particles and one with smaller particles, at different 
radial positions in the protoplanetary disk in contrast to a single 
continuous population from small to large sizes. Our need for 
two different dust grain- size distributions to properly model the 
properties of the disk of the Butterfly Star agrees well with this 
prediction. 

We find that our results are consistent with current models 
for the evolution of dust in circumstellar disks. Our findings of 
grain growth, vertical settling, and radial segregation in the disk 
of the Butterfly Star are common predictions of theoretical mod- 
els. 



8.6. Observability with ALMA 

Figures [17] and [18] show simulated ALMA observations at 
1.3 mm and 894 |im, respectively. They are based on the modeled 
(sub)mm maps of our b est-fit and are created using the CASA 
simulator (Jaeger 2008). Both figures represent the use of the 
full ALMA array with an observing time of 2 hours, in which 
Fig. [TT] reflects configuration 15 and Fig. 1_8 configuration 13. 
The difference between the simulations shown on the left and the 
right in both figures is that for the right one only one single dust 
grain- size distribution with a max = 100 \xm is used whereas the 
simulation on the left exactly represents our best-fit model. The 
dust re-emission in the region that is encased by the 20 % and 
40 % contour level in Fig. 17 left and Fig. _18 left, respectively, 
is strongly dominated by the larger dust grains. Because of the 
much higher angular resolution compared to the SMA and the 
PdBI, the radial extent of this region is cleary visible. The dif- 
ference in the radial brightness structure between the presence of 
only one grain-size distribution or, as we have found in our study, 
two different grain- size distributions in the circumstellar disk of 
the Butterfly Star can be clearly seen in both figures. Thus, ob- 
servations with the full ALMA array will allow us to distinguish 
between the presence of one or two different dust grain- size dis- 
tributions and let us check our current model of the disk of the 
Butterfly Star. Most likely, these observations will allow us to 
give further constraints on the properties of this prominent cir- 
cumstellar disk. 



9. Summary and conclusion 

We have compiled a high-quality data set for the circumstellar 
disk of the Butterfly Star, spanning a wide range of wavelengths. 
We obtained images in the NIR and in the (sub)mm regime, as 
well as photometric data and a spectrum. We have constructed a 
detailed model that allows the interpretation of all of these obser- 
vations with one single set of parameters. A systematic analysis 
of the parameter space allows us to establish strong constraints 
on all the parameters of the model. The conclusions we obtain 
are the following: 

1 . The disk structure is well constrained. The disk extends from 
an inner radius of 0.1 AU up to 300 AU. The scale height of 
the disk is ~ 10 AU at 100 AU and varies with a flaring index 
of -1.14. The surface density exponent is found to be in the 
range of [0.04, 0.48], which is slightly smaller than the range 
found in general for other circumstellar disks (Kitamura 
[eTaLl [20021 1 Andrews & Williams | [20071 jGuilloteau et~aL 
201 lb. The midplane temperature is determined to be 20 K 
atlOOAU. 

2. The dust mass is found to be ~ 9 x 10" 4 M o under the as- 
sumption of spherical grains and Pg ra in = 2.5 gem -3 . The 
dust mass of the region where the large grains are dis- 
tributed is ~ 2.2 x 1O" 4 M and that of the small grains is 
~ 6.8 x 1O" 4 M . With a typical dust-to-gas ratio of 1/100, 
the total disk mass amounts to ~ 9 x 10" 2 M Q , compared to 
a mass of the central star of ~ 1.7 M©. Based on our con- 
straints on the model parameters, the disk of the Butterfly 
Star is found to be non-self-gravitating and is, according to 
the Toomre criterion, gravitationally stable at all radii. 

3. The disk of the Butterfly Star has already entered the first 
phases of planetary formation. We found quantitative evi- 
dence for grain growth, vertical settling, as well as radial 
segregation in this disk. It is not possible to find a coherent 
multi- wavelength model of the disk using a single grain- size 
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Fig. 17. Simulated ALMA observation at 1.3 mm. The contour levels are at 5, 20, 40, 60, and 80 % of the maximum value. 
Left: Simulation based on our best-fit model. Right: Same as left, but using only one grain-size distribution with a m2LX = 100 \xm. 
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Fig. 18. Simulated ALMA observation at 894 |im. The contour levels are at 5, 20, 40, 60, and 80 % of the maximum value. 
Left: Simulation based on our best-fit model. Right: Same as left, but using only one grain-size distribution with a max = 100 \xm. 



4. 



distribution. To do so, the use of at least two different distri- 
butions is necessary. In the outer parts and upper layers of the 
disk, small ISM-sized dust grains can be found. In the inner 
disk regions, larger dust grains with a maximum particle ra- 
dius of up to ~ 100 |im are located. We constrain this region 
where the larger grains are distributed to a radial extent of 
[170, 210] AU and a vertical extent of [4, 7] AU at 100 AU. 
The equivalent scale height of this region amounts to a range 
of [2.5, 4.3] AU and for our best-fit model to 3.7 AU which is 
substantially smaller than the vertical extent and the hydro- 
static scale height (H = 6.9 AU), thus providing strong sup- 
port for vertical settling of the large dust grains in the disk. 
Our findings of dust evolution in this disk are consistent with 
current theoretical models. 

The millimeter spectral index indicates that the dust grains 
in the disk midplane have already grown to sizes larger than 
those found in the ISM, supporting our results. 
Although each individual observation provides valuable in- 
formation on the disk, it is necessary to combine a large 



set of independent observations (SED and spatially re- 
solved images) from different wavelength regimes in a multi- 
wavelength study. This approach allowed us to derive quali- 
tatively new conclusions which were not obvious on the ba- 
sis of individual data sets alone and to strongly reduce model 
degeneracies. 

The large number of observations in different wavelength do- 
mains, as well as our coherent multi- wavelength modeling, make 
the Butterfly Star one of the best studied protoplanetary disks so 
far. With the upcoming completion of the next-generation inter- 
ferometer ALMA, observations using this facility will allow us 
to test our findings of dust evolution, will let us further constrain 
the radial and vertical structure of this fascinating disk, and will 
refine our understanding of the evolution of protoplanetary disks. 
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